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Abstract 

The forming and melting of complementary base 
pairs in RNA duplexes are conformational tran- 
sitions required to accomplish a plethora of bio- 
logical functions. Yet the dynamic steps of these 
transitions have not been quantitatively character- 
ized at the molecular level. In this work, the 
base opening process was first enforced by atom- 
istic pulling simulations and then analyzed with a 
novel reweighting scheme which allowed the free- 
energy profile along any suitable reaction coor- 
dinate, e.g. solvation, to be reconstructed. The 
systematic application of such approach to dif- 
ferent base-pair combinations provides a molecu- 
lar motion picture of helix opening which is val- 
idated by comparison with an extensive set of 
experimental observations and links them to the 
enzyme-dependent unwinding mechanism. The 
RNA intrinsic dynamics disclosed in this work 
could rationalize the directionality observed in 
RNA-processing molecular machineries. 



:> Introduction 



The ability of ribonucleic acid (RNA) to adopt pe- 
culiar three-dimensional structures that mediate a 
variety of biological functions makes it the most 
versatile regulatory factor in the cell.!^ Virtually 
involved in all cellular processing of the genetic 
information, the RNA is able to achieve such a 
functional diversity by adaptively acquiring very 
distinct conformations in response to specific con- 
ditions of the cellular environment.'^' Among the 
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Structural rearrangements engaged by RNA, the 
opening of complementary base pairs is an ubiqui- 
tous process required to accomplish a wide range 
of metabolic activities such as transcription, pre- 
mRNA splicing, ribosome biogenesis or transla- 
tion initiation. ^ In the cell this is usually catalyzed 
by enzymes called RNA helicases which have 
been shaped by the evolution to unwind double- 
stranded (ds) RNA according to its intrinsic dy- 
namic properties."* ^ 

From a molecular standpoint, the opening and 
forming of individual base pairs are fundamental, 
yet poorly understood, events which provide the 
structural framework to large-scale RNA confor- 
mational transitions and folding.^ ^'^^ In this re- 
spect and related to the work presented herein, 
insightful investigations have been reported only 
for short deoxyribonucleic acids (DNAs) in the B- 
form helical geometry.' ^^ Using transition-path 
sampling, Hagan et al^ have fully characterized 
the energetics of (un)pairing for a 5 '-end cytosine. 
However, the mechanism underlying the complete 
opening of the duplex has not been systemati- 
cally faced nor analyzed. Moreover, differences in 
topology and thermodynamic parameters between 
B-form DNA and A-form RNA suggest that the 
mechanism of duplex separation might obey dif- 
ferent rules. 

Recently, combining thermodynamic informa- 
tion with the relative population of unpaired termi- 
nal nucleotides (dangling ends) observed in large 
ribosomal RNA (rRNA) crystal structures, Mo- 
han et al^SI have proposed that stacking and pair- 
ing reactions are not simultaneous, and that 3'- 
single- strand stack leads the base pairing of the 
5 '-strand. Nevertheless, collecting an unbiased 



data set of dangling-end population is not trivial 
and, when viewed in the context of the full ribo- 
somal assemblies, the single stranded (ss) regions 
are seen to interact extensively with other RNA el- 
ements .^51 On top of that, since the formation and 
opening of base pairs is a dynamic process, both 
the ensemble-averaged thermodynamic properties 
and the detailed but static X-ray picture have to be 
complemented with other methods able to directly 
and quantitatively capture the dynamics of the in- 
vestigated event. Likely, this gap will be efficiently 
bridged by ad h oc designed spectroscopic ap- 
proaches. 12121121111 poj- instance, femtosecond time- 
resolved fluorescence spectroscopy is emerging as 
a powerful technique for the quantitative analy- 
sis of base- stacking pattern and base motion,^^ al- 
though its applications to probe RNA dynamics 
are still in their infancy and the method presents 
several limitations.^^ As a matter of fact, the in- 
tegration of spectroscopic approaches with other 
powerful techniques is presently needed to gain 
molecular details on the RNA intrinsic dynam- 
ics. Among the possible methodological choices, 
atomistic simulations^^ allow any base sequence 
to be characterized and all the microscopic param- 
eters to be controlled. Additionally, when com- 
bined with state-of-the-art free-energy methods, 
they can provide an unparalleled perspective on 
the mechanism and dynamics of the biomolecu- 
lar process of interest. As far as the reconstruction 
of free-energy profiles is concerned, the capabil- 
ity of estimating those profiles along any suitable 
reaction coordinate, without any further computa- 
tional cost, would offer researchers a powerful and 
versatile tool enabling both the disclosure of in- 
termediate states and the multifaceted analysis of 
complex conformational transitions. 

With this spirit, here we report an in silico study 
elucidating the mechanism for strand separation 
in the RNA double helix. In particular, we used 
atomistic steered molecular dynamics (MD) simu- 
lations'^ to enforce the unbinding of nucleobases 
into the surrounding explicit water. To allow a sys- 
tematic analysis of different base sequences we de- 
vised a novel Jarzynski-equation-based reweight- 
ing scheme which allowed the free-energy land- 
scape to be reconstructed as a function of different 
reaction coordinates and the unbinding energies 
to be straightforwardly estimated. The computed 



free-energy differences are consistent with exper- 
imental observations and suggest that the strand 
separation mechanism occurs by a stepwise pro- 
cess in which the probability of unbinding of the 
base at the 5' terminus is significantly higher than 
that at the 3' terminus. The biological implica- 
tions of these findings are discussed and related to 
the unwinding mechanism catalyzed by RNA pro- 
cessing machineries. Given the general nature of 
our approach, the introduced methodology can be 
directly applied to analyze a broad range of molec- 
ular unbinding processes. 

Methods 

Throughout the manuscript the following nomen- 
clature will be consistently used to define each el- 
ementary step involving one single base and oc- 
curring during the opening of a closed base pair: 
unpairing is used to define the process undergone 
by a single base for which both Watson-Crick hy- 
drogen bonds and stacking interactions with adja- 
cent bases are broken; unstacking is the process 
breaking the stacking interactions between a dan- 
gling terminal nucleotide and its adjacent bases. 
The opening of a base pair is thus composed by an 
initial unpairing followed by an unstacking ([!]). In 
the manuscript we also use the term unbinding re- 
ferring to both unpairing and unstacking processes 
by no means of specificity. Finally, the strand with 
the 5' terminal (or 3' terminal) nucleobase being 
pulled is referred to as the 5 '-strand (or 3 '-strand). 



System set-up 

We simulated the unpairing and unstacking of 
nucleobases at both 3'- and 5 '-termini in dsR- 
NAs of seauence 5'-ccgggc-3' ^ 5'-ggcccg-3- 

iN/\j, Oi sequence 3'.gGCCCG-5' '^"" 3'-CCGGGC-5' 

([2]). Two sets of data can be obtained from each 
dsRNA thus resulting in four systems with dif- 
ferent combinations of Watson-Crick base pair- 
ing and stacking (|2]\). Both terminal and non- 
terminal base pairs (i.e. a base pair at the ss-ds 
RNA junction) were investigated. Non-terminal 
base pairs showed the same trend in relative sta- 
bility observed for terminal ones, and are reported 
in the Supporting Information (SI). The A-form 



5'-dangling 




Figure 1: Elementary steps involved in the open- 
ing of a base pair. The thermodynamics cycle was 
used to characterize different base pair combina- 
tions. 

dsRNA was built using ASSEMBLE^i and then 
solvated with ~3600 water molecules, 20 Na+ and 
10 Cl^ ions, resulting in an excess salt concentra- 
tion of about 0.15 M. The mobility of added ions 
was fairly diffusive during the simulations. After 
minimization and thermalization, each system (or 
intermediate) was then evolved for 30 ns i n the 
isothermal-isobaric ensemble (300K, lAtm)'2ZE2] 
using the Amber99 force field'^ and TIP3P wa- 
|-gj.i25j Preliminary calculations carried out using 
the recent refinement of the Amber99 force field 
(parmbscO)'2S have shown quantitatively similar 
results in the reconstructed free-energy profiles. 
This is probably due to the poor involvement of 
the refined a and 7 dihedrals during the unbind- 
ing trajectories. Long-range electrostatic interac- 
tions were calculated with the particle mesh Ewald 
method-I^Zl Plain MD and biased steered MD tra- 
jectories were generated with GROMACS 4.0.71^ 
combined with PLUMED 1 2!^ 



Pulling simulations 

The starting configurations for the pulling simu- 
lations were randomly sampled from the corre- 
sponding 30 ns-long runs. The distance between 
the center of mass of two stacked bases ([2^) 
was used as pulled collective variable (CV), and 
thus harmonically restrained to a constant- velocity 
moving point, starting at a position equal to the 
equilibrium average of the CV and pulling it by 




Figure 2: RNA double helix. A) Schematic view 
of the combinations (red dotted boxes) of gua- 
nine (orange) and cytosine (blue) base-pairing and 
-stacking investigated. B) Structural representa- 
tion of the RNA duplex in water; the distance be- 
tween the center of mass (green arrow) of the six- 
membered ring atoms (thick green sticks) of two 
stacked bases was used as collective variable for 
the pulling simulations. 

0.75 nm in 1 .5 ns. This resulted in a biasing poten- 
tial equal to %ag(<?,0 = \[s{q{t))- {sQ + vt)f 
where k = 1200kcal/mol/nm^ is the spring con- 
stant of the restraint, q are the microscopic coordi- 
nates, s{q) is the CV value for those coordinates, 
So is the initial position of the restraint, v is the 
pulling velocity and t is the time. 

The mechanical work done during the process 
was obtained by integrating the force exerted 
on the system along the biased reaction coordi- 
nate. After collecting about 400 realizations for 
each nucleobase-unbinding process, the Jarzynski 
nonequilibrium work theorem ^*^ was exploited to 
discount the dissipated work and to reconstruct the 
free-energy profile as a function of the restraint 
distance (sQ + vt). Although employing Jarzynski's 
equality in principle allows unbiased free-energy 
differences to be estimated, its direct application is 
limited by the number of collectable realizations as 
well as by the complexity of the system-ISIl a typi- 
cal free-energy profile is shown in[3j^ as a function 
of the restraint distance. The blue plot shows how, 
after a steep rise, a series of alternating shoulders 
and local plateaus gradually brought the system to 
higher free-energy states. Moreover, between dis- 
tances ranging from 1 to 1 .2 nm (for the exempli- 
fied system), the profile was strongly dominated 
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Figure 3: Typical nucleobase unbinding process, 
obtained by pulling along the distance between the 
center of mass of two stacked bases. A) Mechan- 
ical work (W) performed (gray plot) and its expo- 
nential average (blue plot) as in Jarzynski equal- 
ity, plotted as a function of the restraint position. 
The distance is practical for biasing the system but 
hardly allowed defining the bound and unbound 
states. B) Number of water molecules coordinat- 
ing the unbinding nucleobase (//wat) as a function 
of time (main panel), and its probability distribu- 
tion (right panel). The coordination with water is 
an useful metrics for identifying the bound and un- 
bound states. 



by an outlier low-work realization and the differ- 
ence in reconstructing the free-energy profile with 
or without the outlier was more than 3kcal/mol 
(see SI). 

Within this framework, there was no clean way 
to automatically detect when the nucleobase had 
reached the unbound configuration, and it was dif- 
ficult to avoid systematic errors in the comparison 
of many profiles with small differences as we were 
interested in. To tackle this problem, we decided 
to analyze our simulations in terms of CVs differ- 



ent from the one used for the pulling. This a pos- 
teriori analysis could be done quickly, as a post- 
processing, and allowed us to choose optimal CVs 
capable of describing in an user-independent man- 
ner all the unbinding events. 

Reweighting scheme 

To project the free-energy landscape on putative 
CVs we devised a proper reweighting scheme. 
Whereas suitable schemes have been proposed to 
reweight other types of nonequilibrium simula- 
tions (e.g. Ref.^, a reweighting algorithm for 
steered MD has not been reported. For a different 
purpose. Hummer and Szabo developed a method 
which enables the reconstruction of the free e nergy 
as a function of the pulled coordinate.'22IHI Here, 
we generalize this scheme so as to compute the 
free energy as a function of any a posteriori cho- 
sen variable. 

Two different sorts of bias affect the steered MD 
trajectories and needed to be removed: (a) the 
nonequilibrium nature of the pulling and (b) the 
presence of artificial harmonic restraints on the 
pulled CV The nonequilibrium bias is removed by 
noticing that the equilibrium probability Peq{q, t), 
for a restraint statically kept in its position at time 
t, can be obtained from the non-equilibrium one 

^neq(<?,0 ^^ observed in the i-th trajectory exploit- 
ing a relation first reported by Crooks: ^^ ^^ 



/'eq(^,0=Ie-/^t^'«-^^^^l/'S(^,0, (1) 

i 

where Wi{t) is the work done on the i-th trajectory 
up to time t and /3 = l/hgT is the inverse thermal 
energy. Here the free energy F{t) represents the 
normalization factor corresponding to the instan- 
taneous position of the moving restraint at time t. 
The bias of the harmonic restraint can then be re- 
moved by applying the weighted-histogram anal- 
ysis method.'23 Whereas weighted histograms are 
traditionally used to combine independent simu- 
lations performed with different static biasing po- 
tentials, here we used it to combine snapshots ob- 
tained at different stages of the pulling, thus writ- 



ing the unbiased equilibrium probability as 



Puiq) 



J^dtP^qiq, 



/o d^^ 



-p[V{q,t)-F{t)] 



(2) 



where T is the length of each pulling simulation. 
Finally, the free energy as a function of an arbi- 
trary, a posteriori chosen CV 5 is defined as F{s) = 
— ksT log J dqPu{q) 5 {s — s{q)). The scheme de- 
scribed so far closely resembles the one used by 
Hummer and Szabo.^^ However, it is conceptu- 
ally different, as here the free energy can be re- 
constructed also with respect to a variable differ- 
ent from the pulled one. Thus, it potentially en- 
ables the disclosure and characterization of other- 
wise hidden features of the investigated process. 
To further simplify the data manipulation and to 
avoid building multidimensional histograms, with 
a further dependence on technical choices such as 
binning size, we recast our approach assigning a 
weight to each of the s amp led configurations, in 
the same spirit as in Ref.^^. After simple manipu- 
lation, the weight can be shown to be equal to 



Wi{t) 



-li[Wi{t)-F{t)] 



lodt'e 



-p[V{qiit),t')-F{t' 



(3) 



The normalization factor for each time, F{t), is 
then computed iteratively up to convergence as 
e-mt) = £./Q^rf?'w,(r>-/5^(?'(0.0. Usually a 

few tens of iterations are enough to converge. 

In summary, in our reweighting scheme we first 
compute the weight of each of the configurations 
saved along the MD simulations from[3j then esti- 
mate free energies as a function of any, a posteriori 
chosen CV as 

F{s) = -kBT\ogY,Mt)Si^- KQ.it))- (4) 



Results 

Using the reweighting scheme outlined in the pre- 
vious Section we were able to investigate several 
order parameters. Since solvent interactions are 
known to affect the conformational state of nucleic 
acids, -^^ we considered the solvation of an unbind- 
ing base as an effective metrics for the progres- 
sion of the underlying process. This choice al- 



lowed defining the unbinding in a manner which 
was totally independent from both the terminus 
and the specific base, and, in our explicit-solvent 
simulation, could be computed as the coordination 
among heavy atoms of the base and water oxygens 
([3p). In this metrics, the bound and unbound states 
could be clearly and unambiguously identified and 
corresponded to approximately harmonic basins. 
Sample free energies computed as a function of the 
number of water molecules coordinating the un- 
pairing base are shown in |4} The free-energy pro- 
file reconstructed along such a reaction coordinate 
is in no way biased by the absolute number of co- 
ordinated water molecules which is merely used to 
distinguish one configuration from the other and to 
properly collect the corresponding weights along 
the simulation as in |4} Then, to compute accu- 
rately the bound/unbound free-energy differences, 
we fit the free-energy profiles with the combina- 
tion of two quadratic f unctions, S^E] 



,-PFi^) 



aj ^e 



lai 



-pFi 
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-1 



-^-/^^^ 



(5) 

where Fi and F2 are the free energies of bound 
and unbound states. Both, when the two states 
were clearly resolvable (e.g. |4| left panels), and 
when the corresponding CV population was more 
overlapped, (e.g.|4[ right panels) the fitting proce- 
dure resulted robust and poorly sensitive to outlier 
work realizations, thus enhancing convergence of 
the results (e.g. the difference in performing the 
fit with or without the outlier low-work realization 
in |3]^ was less than 0.3 kcal/mol, see SI). Fur- 
thermore, this approach showed very stable out- 
comes with respect to the choice of the details in 
the definition of the solvation order parameter, and 
allowed comparing systematically several similar 
situations without incurring of large statistical er- 
rors or, worst, human biases in the interpretation 
of the results. 

Having an optimized statistical-mechanics tool 
able to provide free-energy differences in a flexible 
and automatic manner, we pursued a systematic 
step-by- step approach to investigate the feasibil- 
ity of different opening paths for the four possible 
combinations of G-C base stacking and Watson- 
Crick pairing. The general procedure, as outlined 
in [T| relied on two subsequent steps: first, the 




Figure 4: Reconstruction of the free-energy pro- 
file as a function of the number of water molecules 
(A'^wat) surrounding the unbinding base. In the left 
and right panels are shown typical free energy pro- 
files (red dots with error bar) for the unpairing of 
a 3'-strand guanine and 5'-strand cytosine, respec- 
tively. The quadratic potentials obtained from the 
double-well fitting are shown in light blue color, 
whereas their combination [|5| is in blue. Under- 
neath each panel, the unnormalized population of 
the CV is also shown. 

Watson-Crick base pair was partially opened by 
the unpairing of the base on either the 5' or 3' 
terminus; second, the resulting dangling interme- 
diate, on the 3' or 5' terminus respectively, was 
unstacked and the base pair opening completed. 

Table 1: Context-dependent base-unbinding free 
energy (kcal/mol) corresponding to the elementary 
steps shown infTl 

Construct Opening steps 
n 3b 5b 5a 3a 



(1) 



5'-CG.. 
'3'-GC.. 

'^^^3'-GG.. 
/on5'-GG.. 
v-^''3'-CC.. 



7.6 0.9 2.6 5.6 

8.9 0.3 3.9 4.5 

7.0 2.2 4.9 3.9 

7.4 2.0 5.7 2.8 



The relative stability of putative intermediates 
involved in the opening of a base pair was esti- 
mated from the individual base-unbinding free en- 
ergies ([T]and[T]). For all the considered combina- 
tions, the difference in basepair-opening free en- 
ergy computed biasing the system along path (a) 
and (b) in[T]was lower than Ikcal/mol. For sake of 
clarity, it should be reminded that a finite number 



of unidirectional pulling simulations performed 
within a Jarzynski-like scheme are known to pro- 
vide overestimates of absolute free-energy differ- 
ences.^^ However, highly accurate estimates of 
unbinding constants were not needed to character- 
ize the strand separation mechanism and the free- 
energy differences we estimated were exploited as 
a quantitative tool to assess the relative stability of 
different configurations. 



Comparison with experiments 

Below we discuss the results of the first and sec- 
ond unbinding steps ([T]), and compare t hem w ith 
crystal structures conformer distributions, T^EII rel- 
ative population of stacked/unstacked bases de- 
tected by femtosecond time -resolved fluorescence 
spectroscopy,S3 and thermody namic data based 
on dsRNA melting experiments.^ ^ | 44 | 45 ] j^^^ consis- 
tency with experimental observations and the ca- 
pability of our simulations to complement those 
results are highlighted. 

The more general outcome arising from the com- 
parison of free-energy differences is that the paired 
base on the 5' terminus always interacted more 
weakly than the complementary base on the 3 ' ter- 
minus (steps 5a, 3b in [T] and [T]). This could be 
directly related to the A-form helical geometry of 
RNA in which the bases at the 5' end of a ss-ds 
junction are less buried into the neighboring envi- 
ronment and expose a wider portion of their sur- 
face to water molecules, thus facilitating fraying 
events. The different stability of the nucleobase 
on the 5' terminus can be reflected in the prob- 
ability of observing a certain type of blunt clos- 
ing base pair at ss-ds junctions. In this context, 
the stronger interaction was estimated for the 5'- 
guanine in ^/^q (construct 4) which was ~1.7 
kcal/mol weaker than the complementary cytosine 
on the 3'-terminus. Accordingly, the combina- 
tion j/j^q" (construct 4) is the most abundant clos- 
ing base-pair pattern observed at ss-ds junctions 
in large rRNA crystal structures .l^^ It can be fur- 
ther noticed that among the dangling ends (steps 
5b and 3a in [T] and [T]) the most stable ones are 
those on the 3' terminus, consistently with ul- 
trafast spectroscopy experiments which have de- 
tected a large subpopulation of stacked conform- 



ers for a 3 '-dangling fluorescent purine probe, 
while only a relatively small one for a 5 '-dangling 
purine probe.'^ In particular, we found the most 

stable 3'-dangling end in construct 1 i^'.Qc)^ 
which has also been counted as the most com- 
mon dangling end pattern in rRNA crystal struc- 
tures. ^^ Further agreement can be found consid- 
ering dsRNA optical melting experiments which 
have shown that single-nucleotides overhanging at 
3 '-ends of an RNA helix increase the stability of 
the duplex in a sequence-dependent manner. No- 
tably, such a stabilization has been interpreted as 
the capability of the 3' dangling ends to stack 
over the hydrogen bonds of the closing base pair 
protecting them from water exchange. ^^ Those 
3 '-dangling bases which are more likely stacked 
would thus provide a larger contribution to duplex 
stabilization. In this light, the trend that we ob- 
served in the unstacking energies of the four dan- 
gling constructs (3':2:>3'-Ga >3'S::>3'-cg;.) is 
in agreement with duplex stabilization observed 
in dsRNA melting experiments.^ ^ |44 | 45 ] j^. gj^guld be 
remarked that the duplex stabilization induced by 
5 '-dangling ends might not reflect the stacking en- 
ergy of the dangling end itself because of its small 
overlap with the hydrogen bonds of the closing 
base pair. Further discussion on the comparison 
of computed and experimental dangling-end sta- 
bilities can be found in the SI. 

Summarizing, the unbinding of the base on the 
5 '-strand was, in all the considered cases, favored 
over the unbinding of the complementary 3 '-strand 
base. Whereas the relative probability of 3'- and 
5 '-unbinding event can be modulated by the se- 
quence, the general trend remains unchanged. 

From a structural standpoint, the unpairing could 
proceed through two qualitatively different paths: 
one in which the twisting and breaking of Watson- 
Crick hydrogen bonds occurred before the rupture 
of stacking interactions; the other, in which the 
unbinding followed the concerted rupture of both 
hydrogen bonds and stacking interactions with, in 
some cases, the unbinding base stacking over the 
dangling end of the opposite strand ([5]). Similar 
unbinding geometries have also been described in 
other studies.13E3 in our simulations, these inter- 
mediate states occurred with a context-dependent 
frequency. For instance, along the unbinding path- 



way the 5 '-terminal guanine in y'_QQ" (construct 
4) had ~15% of probability to stack upon the 
3 '-dangling cytosine of the opposite strand. In- 
terestingly, this probability dropped to ~3% in 
3," (^(-.'(construct 3). Such an inter-strand stacking 
pattern (|5} panel t2), exchanging with the conven- 
tional stacking of a paired base, could account for 
the delayed quenching of fluorescence detected by 
ultrafast fluorescence spectroscopy for a construct 



similar to 



5'-GC.. 
3'-CG.. 



(construct 4). 
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Figure 5: Snapshots sampled from the opening of 
a base pair. The unbinding base (here a 5 '-terminal 
guanine in construct 4) could transiently stack over 
the dangling end of the opposite strand (12). 

As final experimental evidences which corrob- 
orate our results, Xia and co-workers ^^ have re- 
ported that the dynamic behavior of a 3 '-terminal 
purine is not affected by the presence of the op- 
posite complementary base. Viceversa, the con- 
formational dynamics of a 5 '-terminal purine is 
drastically influenced by the presence of an oppo- 
site 3 '-terminal pyrimidine which would be likely 
stacked and potentially able to shift the population 
of the complementary 5 '-terminal base towards a 
paired and stacked ensemble. Consistently with 
our systematic study, these data depict the forma- 
tion of a stable base pair as generally driven by 
the stacking of the 3 '-terminal base, and then by 
the energy gained by the system from both the 
stacking of the 5 '-terminal base and Watson-Crick 
hydrogen-bonds. 



tracking tunneu 



Figure 6: Model for RNA unwinding catalyzed by NS3 helicase. A) The single-stranded 3' terminus is 
loaded into the tracking tunnel. B) The 5' terminus is mechanically displaced by the helix opener. C) The 
helicase proceeds with 3' to 5' directionality, displacing the 5 '-strand. 



Biological implications 

The motion picture of duplex separation emerg- 
ing from the outcome of our simulations comple- 
ments and augments with dynamic details and en- 
ergetic considerations the helix propagation model 
based on the analysis of static 3D structures.'^ Our 
computations link experimental data from differ- 
ent fields creating a common reading frame among 
them. Taken together, these results suggest that 
RNA unwinding occurs by a stepwise process in 
which the probability of unbinding of the base on 
the 5' strand is significantly higher than that on 
the 3' strand. What could be the biological impli- 
cations of this finding? 

When considering the RNA as the substrate of 
molecular motors such as helicases and other re- 
modeling enzymes, the results could likely be in- 
terpreted from an evolutionary point of view which 
could allow deciphering the basis of the evolution- 
ary pressure responsible for the unwinding mech- 
anism catalyzed by RNA-duplex processing en- 
zymes. 

The RNA unwinding catalyzed by helicases is 
coupled to adenosine triphosphate (ATP) binding 
and hydrolysis. The underlying mechanism would 
reasonably minimize the use of ATP, especially in 
a low-nutrient environment. Provided that the in- 
trinsic RNA dynamics implies that at ss-ds junc- 
tions the unbinding of the 5 '-strand base is fa- 
vored over the unbinding of the complementary 
3 '-strand base, an enzymatic unwinding model 
would include a mechanism in which the separa- 
tion of the two complementary strands is accom- 



plished by acting on the weakest portion, i.e. the 
5 '-strand base. Thus, an ancestral enzyme using 
the 3 '-strand as running track rail (with 3' to 5'- 
directionality) without perturbing its conformation 
and causing the displacement of the 5 '-strand by 
mechanical exclusion could satisfy some energy- 
saving requirements. 

The viral RNA helicase NS3 of hepatitis C virus, 
which is a prototypical DEx(H/D) RNA helicases 
essential for viral replication, could satisfy the 
above mentioned requirements."^ NS3 is a po- 
tentially relevant drug target and has been struc- 
turally and functionally characterized in various 
contexts .'^^"511 n unwinds duplexes by first load- 
ing onto a single-stranded 3 '-terminus region and 
then processively translocating with 3' to 5' direc- 
tionality along this loading strand, thereby peeling 
off the complementary 5 '-strand bases. In particu- 
lar, the 3 '-strand would migrate through a tracking 
ssRNA tunnel running within the protein whereas 
the complementary 5 '-strand is forced towards the 
back o f the protein by the "helix opener" hair- 
pin (|6|) .122123 In light of the free-energy calcula- 
tions discussed above, it could be suggested that 
this mechanism has been optimized according to 
the intrinsic RNA unwinding dynamics disclosed 
in this work. 

Arguably, the processing machineries are being 
constantly shaped by the evolutionary pressure of a 
plethora of (often unknown) factors contributing to 
the optimization of metabolism in the whole living 
system, rather than to the local biochemical pro- 
cess. As a consequence, the preference for a well- 
defined RNA processing directionality cannot be 



ubiquitously observed.^El 

We speculate that all the biochemical processes 
involving RNA in which directionality plays a role 
(e.g. transcription), could be related to the ener- 
getics of RNA double hehx forming and fraying 
discussed in this Article. 

Conclusions 

This study lays down the basis for the molecular- 
level understanding of intrinsic RNA dynamics 
and its role in function. The asymmetric behav- 
ior of the 3'- and 5 '-strand could be responsible 
for the directionality observed in RNA processing. 
From a computational perspective, the approach 
we introduced can be generalized to analyze any 
kind of (un)binding event. Indeed, it allowed the 
free-energy landscape to be reconstructed along 
different reaction coordinates and the unbinding 
energies to be easily computed in an automatic and 
user-independent manner, therefore removing sta- 
tistical and human biases. We foresee the applica- 
tion of our approach to a wider range of molecular 
systems, including the typical ligand-target com- 
plex faced in drug discovery .I^S 
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